Il dataset che andremo ad analizzare contiene le osservazioni giornaliere del meteo, ricavate dalle rilevazioni effettuate da varie stazioni metereologiche Australiane. L'obiettivo di questo progetto è di addestrare un classificatore binario, in modo che riesca a predire se il giorno dopo pioverà oppure no.
Le analisi del nostro progetto si dividono, principalmente, in tre fasi:
Le prime due sono le fasi principali, nelle quali creiamo la versione finale del dataset e tramite esso otteniamo vari modelli di classificazione. L'ultima fase, è utile per modificare degli aspetti del dataset in modo da ottenere dei risultati di classificazione migliori. Il post-processing viene effettuato solamente sui classificatori che hanno ottenuto i risultati migliori. Quindi, questa fase è come se fosse una sotto-sezione della fase di classificazione.
1.1 Informazioni generali del dataset
1.2 Analisi, statistiche e correlazione tra gli attributi
1.3 Discretizzazione
1.4 Dataset finale
2.1 Clustering (Extra)
2.2 Vari modelli di classificazione
2.2.1 Training e Performance dei modelli
2.2.2 Ensemble
2.2.3 Post-Processing
2.3 Rete neurale
2.3.1 Training e Performance
2.3.2 Post-Processing
In questa prima sezione andremo ad analizzare il dataset e ad effettuare delle statistiche sugli attributi. Cercando di capire quali attributi sono importanti, quali bisogna modificare, affinchè si possano ottenere dei buoni modelli di classificazione, e quali di essi vadano eliminati.
Per iniziare, importo il dataset e cerco di ottenere semplici informazioni riguardo la dimensione e statistiche sulle colonne.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('seaborn-notebook')
df = pd.read_csv('weatherAUS.csv')
df.head(5)
In totale sono presenti 24 attributi, la cui descrizione di ognuno di essi è riportata qui sotto:
Ora controllo la dimensione del dataset, per capire con quante entry abbiamo a che fare.
df.shape
Verifico quante entries sono presenti per ciascun attributo, così da verificare la presenza di valori nulli.
df.info()
Descrizione generale delle statistiche del dataset.
df.describe()
Le statistiche di risk_MM sono quasi equivalenti a Rainfall, meglio verificare una loro possibile correlazione. Inoltre il 75% dei dati di rainfall è prossimo allo 0, potrebbe essere utile discretizzarlo. Prima di fare tutto questo, controllo le colonne con il maggior numero di valori nulli. Ovvero, verifico in percentuale quanti nan sono presenti nelle corrispondenti colonne.
pd.set_option("max_columns",24)
pd.DataFrame(df.isnull().sum()/df.shape[0]).T
Evaporation e sunshine presentano troppi nan, meglio eliminarli per evitare di falsare le future analisi. Droperò anche le righe che presentano tutti valori null.
df.dropna(axis = 0, how = 'all')
df = df.drop('Evaporation', axis = 1)
df = df.drop('Sunshine', axis = 1)
df.shape
Le prime analisi, verranno effettuate sugli attributi che presentano maggiori valori nun e tra quelli che presentavano statistiche simili. Infatti verifico la correlazione tra Rainfall e tutti gli altri attributi, ma in particolare verifico, se è presente una forte correlazione con Risk_MM. Dato che in precedenza questi due attributi presentavano le stesse statistiche.
mtxcorr = df.corr()
mtxcorr["Rainfall"].sort_values(ascending=False)
Correlazione bassa, però i valori differiscono di poco; considerando le statistiche quasi equivalenti (cambia la media ma di poco). Quindi ho deciso di dropparla, considerando che anche la descrizione di Kaggle consiglia di farlo; essendo una collona aggiunta quando hanno messo insieme i dati, e non è chiaro come l'abbiano ricavata. Meglio non rischiare ed eliminiamola.
df=df.drop('RISK_MM',axis=1)
df.shape
Altri due attributi presentano un'elevata percentuale di valori nulli, 40% circa ,ma sono solamente cloud3pm e cloud9am. Leggendo la descrizione su kaggle il valore delle nuvole è rappresentato in scale okta, nella quale nan indica nil significant cloud, ovvero che ci sono le nuvole, ma la loro presenza non è significativa. Quindi nel nostro caso, una presenza non significativa di nuvole assume lo stesso significato di 0 (niente nuvole), nel nostro caso della classificazione. Ipotizziamo che nan non sia un campionamento non effetuato, ma un valore della scala okta come detto poc'anzi. Se per caso la nostra ipotesi si dovesse rivelare errata, ovvero il valore della classificazione dovesse risultare non soddisfacente, provvederemo ad eliminare questi attributi e verificheremo se la classificazione avrà dei miglioramenti oppure no.
df.Cloud3pm.unique()
df.Cloud9am.unique()
I valori assunti dalle due colonne rispecchiano i valori della scala okta. Prima di rimpiazzare i nan con 0, vorrei analizzare i casi in cui è presente il numero 9. Perchè il 9 indica cielo oscurato da nebbia o altri fenomeni atmosferici. Quindi questo presenta delle anomalie atmosferiche, che potrebbero essere delle anomalie anche nel nostro dataset o dei semplici utliers.
df[df['Cloud3pm']==9]
df[df['Cloud9am']==9]
Potrebbero essere degli outliers però il valore 9 è presente solo a uno dei due degli attributi e mai in contemporanea, quindi probabilmente non falsificano le analisi.C'è addirittura una correlazione, che quando è presente il 9 è sempre NO l'etichetta della piogga. Anche se non è significativa come correlazione dato che sono solo 3 casi.
#riempio i nan con 0
df['Cloud3pm'].fillna(0,inplace=True)
df['Cloud9am'].fillna(0,inplace=True)
df.info()
Le due colonne riguardanti le nuvole non presentano più valori nulli. 142193 valori non-null come le entries del dataset.
Adesso effettuo delle analisi statistiche, riguradanti la probabilità di pioggia quando è presente un determinato valore degli attributi Cloud9am e Cloud3pm.
# vediamo che succede quando c'è presenza di nuvole
df_Cloud = df.groupby(['Cloud9am', 'RainTomorrow'])[['Date']].count()
df_Cloud = df_Cloud.rename({'Date': 'count'}, axis=1)
df_Cloud['prob'] = df_Cloud['count'] / df.groupby('Cloud9am').count()['RainTomorrow']
df_Cloud
Non è molto correlato con la pioggia, l'unico dato significativo l'abbiamo con 8, che indica pioggia quasi al 50%; già con 7 piove al 32% e man mano che diminuisce il valore diminuisce sempre di più la probabilità di pioggia.
df_Cloud3 = df.groupby(['Cloud3pm', 'RainTomorrow'])[['Date']].count()
df_Cloud3 = df_Cloud3.rename({'Date': 'count'}, axis=1)
df_Cloud3['prob'] = df_Cloud3['count'] / df.groupby('Cloud3pm').count()['RainTomorrow']
df_Cloud3
Statistiche simili rispetto a prima anche con l'attributo Cloud3pm, anche se la probabilità di pioggià è leggermente più alta.
Dato che i valori delle nuvole sono correlati tra di loro, presentano le stesse statistiche e la stessa probabilità di pioggia, faccio un'unica colonna che presenta la media dei due valori. Sarà dato maggiore credito alle nuvole presenti alle 9pm.
df['Cloud3pm']= ((df.Cloud3pm*8)+(df.Cloud9am*2))//10 #calcolo pesato 80% di importanza alle 3 del pomeriggio rispetto alle 8 del mattino
df['Cloud3pm']
df=df.rename(columns={'Cloud3pm':'Cloud'})
df=df.drop('Cloud9am',axis=1)
df.head(5)
Sono riuscito a mettere solo una colonna. Ora verifico se la correlazione è rimasta equivalente.
# vediamo che succede quando c'è presenza di nuvole
df_Cloud = df.groupby(['Cloud', 'RainTomorrow'])[['Date']].count()
df_Cloud = df_Cloud.rename({'Date': 'count'}, axis=1)
df_Cloud['prob'] = df_Cloud['count'] / df.groupby('Cloud').count()['RainTomorrow']
df_Cloud
Correlazione rimasta equivalente.
Ora mi occupo di un'altra colonna che presenta un basso numero di valori nulli, ma potrebbe rivelarsi un attributo fondamentale. Sto parlando della colonna RainToday, alla quale è possibile sopperire alla mancanza di valori, semplicemente andando a inserire nelle celle vuote il valore della colonna di RainTomorrow del giorno precedente. Ovviamente ho controllato che il primo dato, ovvero il primo giorno, non avesse valori nulli, perchè in quel caso sarebbe impossibile sostituire il valore non avendo un giorno precedente da controllare.
#sostituisco i valori null di rain today controllando il valore di RainTomorrow precedente
index=df.loc[df['RainToday'].isnull(),:].index
df['RainToday'][index]=df.loc[index-1,'RainTomorrow']
df.info()
Discretizzo i valori di Rain today e rain tomorrow in 1 e 0, così questi attributi passano da nominali a categorici. Effettuo questo cambio, perchè il calcolatore e gli algoritmi non lavorano bene con gli attributi nominali.
df['RainTomorrow'].replace({'Yes': 1,'No': 0},inplace=True)
df['RainToday'].replace({'Yes': 1,'No': 0},inplace=True)
df.head()
Controlliamo la relazione tra raintoday e rain tomorrow. Se per caso la pioggia del giorno corrente, aumenta la probabilità di pioggia del giorno seguente.
df_today = df.groupby(['RainToday', 'RainTomorrow'])[['Date']].count()
df_today = df_today.rename({'Date': 'count'}, axis=1)
df_today['prob'] = df_today['count'] / df.groupby('RainToday').count()['RainTomorrow']
rain = df_today.loc[1]
no_rain= df_today.loc[0]
fig, axes = plt.subplots(1,2)
axes[0].bar(1,rain.loc[0]['prob'], label='RainTomorrow: no',align='edge',width=0.6,alpha=0.5)
axes[0].bar(1,rain.loc[1]['prob'], label='RainTomorrow: yes', align='edge',width=-0.6,alpha=0.5)
axes[1].bar(1,no_rain.loc[0]['prob'], label='RainTomorrow: no', align='edge',width=0.6,alpha=0.5)
axes[1].bar(1,no_rain.loc[1]['prob'], label='RainTomorrow: yes', align='edge',width=-0.6,alpha=0.5)
axes[0].set_title('RainToday: yes')
axes[1].set_title('RainToday: no')
axes[0].set_ylabel('Probabilità')
axes[1].set_ylabel('Probabilità')
axes[0].set_xticks([])
axes[1].set_xticks([])
axes[0].legend()
axes[1].legend()
fig.set_size_inches(14, 8)
plt.show()
RainToday influisce molto rain tomorrow nel caso in cui non piove, con 80% di probabilità che non piove nemmeno il giorno seguente.Se piove, il giorno seguente avremo una probabilità del 45% che piova. Influisce più di quanto pensassi.
Verifico la correlazione tra TUTTI gli attributi. Così verifico se qualche attributo è la combinazione di altri. Prima plotto i grafici e poi per casi specifici li verifico uno ad uno.
Controllo per primo la matrice di correlazione per tutti gli attributi. Molto probabilmente non si capirà molto, infatti dopo la vedo per singoli attributi.
plt.figure(figsize=(12,8))
sns.heatmap(mtxcorr,annot=True,cmap='bone',linewidths=0.25)
Già da questa matrice si nota la correlazione tra gli attributi che operano sugli stessi settori. Ad esempio notiamo la forte correlazione tra gli attributi che riguardano la temperatura, correlazione del 90% per alcuni di essi. La correlazione tra gli attributi della pressione, 96% di correlazione, e così via. Per una maggiore chiarezza verifico la correlazione degli attributi che rappresentano le stesse grandezze fisiche, tramite delle scatter matrix.
N.B. Non è presente una scatter matrix tra tutte le possibili coppie di attributi, dato che la macchina non riesce a computarla, crasha spesso se provo a farla di tutti.
Per primo controllo le ralzioni tra gli attributi riguardanti la tamperatura.
from pandas.plotting import scatter_matrix
attributes = ['MinTemp','MaxTemp','Temp9am','Temp3pm','RainToday','RainTomorrow']
scatter_matrix(df[attributes], figsize=(36,24),alpha=0.2)
plt.show()
g= sns.pairplot(df,
x_vars = attributes[:-2],
y_vars = attributes[:-2],
hue='RainTomorrow',
diag_kind='kde',
diag_kws={'alpha':0.2},
)
g.fig.set_size_inches(24,12)
Dalla figura si nota una forte relazione tra MaxTemp e Temp3pm, sembra quasi una linea retta. L'altra forte relazione è tra Temp9am e MinTemp, che è quasi evidente come la correlazione precedente. Notiamo che la piovosità del giorno stesso non influisce molto sulla piovosità del giorno dopo, però la tenamo in considerazione, per il momento, perchè potrebbe essere utile o correlata agli altri attributi. Per ora elimino le due colonne aggiuntive della temperatura.
df=df.drop("Temp3pm",axis=1)
df=df.drop("Temp9am",axis=1)
df.shape
Adesso controllo gli attributi adibiti all'umidità.
attributes = ['Humidity9am','Humidity3pm','RainToday','RainTomorrow']
scatter_matrix(df[attributes], figsize=(36,24),alpha=0.2)
plt.show()
g= sns.pairplot(df,
x_vars = attributes[:-2],
y_vars = attributes[:-2],
hue='RainTomorrow',
diag_kind='kde',
diag_kws={'alpha':0.2},
)
g.fig.set_size_inches(24,12)
Forte correlazione tra le due
mtxcorr = df.corr()
mtxcorr["Humidity3pm"].sort_values(ascending=False)
Anche in questo caso decido di unire le due colonne in una. Anche questa volta effettuo una media pesata dando più importanza all'attributo Humidity3pm, in quanto è una rilevazione che temporalmente è più vicina al giorno seguente, rispetto a quella effettuata alle 9 della mattina.
df['Humidity3pm']= (df.Humidity3pm*8+df.Humidity9am*2)/10
df=df.rename(columns={'Humidity3pm':'Humidity'})
df=df.drop('Humidity9am',axis=1)
#controllo se la correlazione con rain tomorrow è rimasta uguale
mtxcorr = df.corr()
mtxcorr["Humidity"].sort_values(ascending=False)
Abbiamo aumentato la correlazione con l'attributo significativo per la classificazione, ovvero RainTomorrow. Quindi l'unine delle due colonne è stata eseguita correttamente.
Adesso verifico la pressione.
attributes = ['Pressure9am','Pressure3pm','RainToday','RainTomorrow']
scatter_matrix(df[attributes], figsize=(36,24),alpha=0.2)
plt.show()
g= sns.pairplot(df,
x_vars = attributes[:-2],
y_vars = attributes[:-2],
hue='RainTomorrow',
diag_kind='kde',
diag_kws={'alpha':0.2},
)
g.fig.set_size_inches(24,12)
Anche in questo caso forte correlazione tra le due colonne.
mtxcorr = df.corr()
mtxcorr["Pressure3pm"].sort_values(ascending=False)
Attributi quasi equivalenti, per evitare questa ridondanza li unisco in una sola colonna effettuando la media aritmetica dato che i due valori sono quasi uguali
df['Pressure3pm']= (df.Pressure3pm+df.Pressure9am)/2
df=df.rename(columns={'Pressure3pm':'Pressure'})
df=df.drop('Pressure9am',axis=1)
df.shape
Ultimo controllo sugli attributi del vento.
attributes = ['WindGustSpeed','WindSpeed9am','WindSpeed3pm','RainToday','RainTomorrow']
scatter_matrix(df[attributes], figsize=(36,24),alpha=0.2)
plt.show()
g= sns.pairplot(df,
x_vars = attributes[:-2],
y_vars = attributes[:-2],
hue='RainTomorrow',
diag_kind='kde',
diag_kws={'alpha':0.2},
)
g.fig.set_size_inches(24,12)
mtxcorr = df.corr()
mtxcorr["WindGustSpeed"].sort_values(ascending=False)
SI nota che entrambi gli attributi WindSpeed sono abbastanza correlati con WindGustSpeed, il che indica che combinati potrebbero dare WindGustSpeed, dato che quest'ultimo rappresenta il cambio di intesità della velocita del vento durante la giornata. Inoltre WindSpeed3pm e WindSpeed9am, presentano una bassa correlazione con rainTomorrow rispetto a WindGustSpeed, il quale è correlato anche alla pressione. Quindi ho deciso di eliminare queste due colonne, comprese anche le colonne relative alle loro direzioni del vento.
df=df.drop('WindSpeed9am',axis=1)
df=df.drop('WindSpeed3pm',axis=1)
df=df.drop('WindDir9am',axis=1)
df=df.drop('WindDir3pm',axis=1)
df.shape
Adesso passo ai dati non categorici, iniziando da data. Inanzitutto controllo fino a dove inizia e fino a dove finisce il periodo di osservazione, controllo anche la presenza di date ripetute.In seguito, la prima analisi statistica che effettuerò riguarda l'andamento generico della piovosità.
df_sort=df.sort_values(by=['Date'],ascending=True)
print("Inizio osservazione: "+str(df_sort['Date'].iloc[0]))
print("Fine osservazione: "+str(df_sort['Date'].iloc[len(df)-1]))
print("Numero giorni (approssimato) di osservazione: "+str((2017-2007)*365-30*5))
I giorni del periodo di osservazione sono effettivamente minori del numero di entries, quindi probabilmente le date sono al loro volta divise per zona. Ovvero a una dato corrispondono più location. Provo a confermare questa ipotesi.
df_check = df_sort.groupby(['Location'])[['Date']].count()
df_check = df_check.rename({'Date': 'NumGiorniOsservati'}, axis=1)
df_check
dateTotali=df_check.NumGiorniOsservati.values.sum()
if dateTotali==df['Date'].count():
print("Ipotesi Confermata!")
else :
print("Ipotesi Errata!")
La nostra ipotesi è stata confermata, perchè il numero di giorni osservati combiacia più o meno con il numero di giorni che passano nel periodo di osservazione. Alcune località però presentano un numero significativamente più basso di valori di data, il che significa che molte zone, come Uluru,Katherine etc. non hanno riportato la loro osservazione per lunghi periodi.
Siccome non c'è uniformita nel numero di osservazioni, proverò a effettuare delle statistiche in base ai mesi in cui ha piovuto di più, in termini di giorni di pioggia, e i mesi in cui è scesa più pioggia in termini di ml.
#Faccio una copia del dataset per suddividere la data in mesi
import re
df_copy=df.copy()
mesi=['Gennaio','Febbraio','Marzo','Aprile','Maggio','Giugno','Luglio','Agosto','Settembre','Ottobre','Novembre','Dicembre']
df_copy.Date.replace({'[0-9]+-01-[0-9]+':0,'[0-9]+-02-[0-9]+':1,'[0-9]+-03-[0-9]+':2,'[0-9]+-04-[0-9]+':3,
'[0-9]+-05-[0-9]+':4,'[0-9]+-06-[0-9]+':5,'[0-9]+-07-[0-9]+':6,'[0-9]+-08-[0-9]+':7,
'[0-9]+-09-[0-9]+':8,'[0-9]+-10-[0-9]+':9,'[0-9]+-11-[0-9]+':10,'[0-9]+-12-[0-9]+':11},
regex=True,inplace=True)
df_mesi=df_copy.groupby(['Date'])[['RainToday']].sum()//(30*10)
df_mesi['MediaPioggiaCaduta']=df_copy.groupby(['Date']).mean()['Rainfall']
fig, ax = plt.subplots(1,1)
fig.set_size_inches(9,9)
ax.plot(df_mesi['RainToday'],'ko--')
ax.set_xticks([i for i in range(12) ])
xlabels = ax.set_xticklabels(mesi,rotation=30, fontsize=12)
ax.set_title('Giorni di pioggia medi per ogni mese')
_ = ax.set_ylabel('Numero Giorni')
ig, ax = plt.subplots(1,1)
fig.set_size_inches(9,9)
ax.plot(df_mesi['MediaPioggiaCaduta'],'ko--')
ax.set_xticks([i for i in range(12) ])
xlabels = ax.set_xticklabels(mesi,rotation=30, fontsize=12)
ax.set_title('Quantità di pioggia media per ogni mese')
_ = ax.set_ylabel('mm di pioggia')
Dalle statistiche si nota che i mesi in cui è piovuto più giorni, coprono il periodo da maggio a luglio. Mentre gli altri mesi è piovuto lo stesso numero di giorni circa, da 2200 a 2500 giorni. Nonostante questo, è caduta più pioggia, in termine di ml, da Gennaio fino a Marzo, fatta eccezione per Giugno. A noi ovviamente interessa il numero di giorni in cui è piovuto e non la quantità. ANche se questi non sono dati molto significativi, considerando che il minor numero di giorni di pioggia medi è 7, mentre il massimo è 11. Quindi i giorni di pioggia sono ben distribuiti nell'arco dei mesi.
Dato che i periodi di pioggia sono divisi per blocchi di mesi, provo a dividerli per stagioni. In modo da raggruppare i mesi in cui piove di iù in un'unica batch.
df_copy=df.copy()
stagioni=['Estate','Autunno','Inverno','Primavera']
df_copy.Date.replace({'[0-9]+-(12|0[1-2])-[0-9]+':0,'[0-9]+-0[3-5]-[0-9]+':1,'[0-9]+-0[6-8]-[0-9]+':2,'[0-9]+-(09|1[0-1])-[0-9]+':3},
regex=True,inplace=True)
df_stagioni=df_copy.groupby(['Date'])[['RainToday']].sum()//(30*3*10)
df_stagioni['MediaPioggiaCaduta']=df_copy.groupby(['Date']).mean()['Rainfall']
fig, ax = plt.subplots(1,1)
fig.set_size_inches(9,9)
ax.plot(df_stagioni['RainToday'],'ko--')
ax.set_xticks([i for i in range(4) ])
xlabels = ax.set_xticklabels(stagioni,rotation=30, fontsize=12)
ax.set_title('Giorni di pioggia medi per ogni stagione')
_ = ax.set_ylabel('Numero Giorni')
ig, ax = plt.subplots(1,1)
fig.set_size_inches(9,9)
ax.plot(df_stagioni['MediaPioggiaCaduta'],'ko--')
ax.set_xticks([i for i in range(4) ])
xlabels = ax.set_xticklabels(stagioni,rotation=30, fontsize=12)
ax.set_title('Quantità di pioggia media per ogni stagione')
_ = ax.set_ylabel('mm di pioggia')
Effettuando la statistica sulle stagioni non abbiamo ottenuto dei dati considerevoli rispetto all'analisi effettuata sui mesi. Anzi adesso ci siamo accorti che non c'è una grossa differenza tra le varie stagioni, i giorni di pioggia medi durante l'anno cambiano di poco. L'unica differenza rigurda l'ammontare di pioggia che cade, d'estate piove di meno in termini di giorni ma di più in termini di mm. Ma tutto ciò non è rilevante per gli obiettivi del nostro tipo di classificazione.
df.Location.unique()
Tutte queste zone andando a controllare su internet, fanno parte della zona del sud est dell'australia. COnsiderando che copre solo quella zona e che in alcune località mancano parecchie date, si è deciso di eliminare le località, perchè potrebbero causare overffitting, basare del training su di esse. Dato che nel momento in cui uscisse una nuova zona potrebbe essere falsata la classificazione, se si basa sull'attributo della location. Inoltre vorrei droppare anche le date, perchè, come visto dai grafici, in termine di piovosità non c'è molta differenza tra le varie stagioni o mesi. Infine, molte di esse presentano dei valori mancanti, quindi alcune statistiche sono falsate.
df=df.drop('Location',axis=1)
df=df.drop('Date',axis=1)
df.shape
Ultimo controllo sull'attributo Wind Gust Direction.
df.WindGustDir.unique()
df_Dir = df.groupby(['WindGustDir', 'RainTomorrow'])[['WindGustSpeed']].count()
df_Dir = df_Dir.rename({'WindGustSpeed': 'count'}, axis=1)
df_Dir['prob'] = df_Dir['count'] / df.groupby('WindGustDir').count()['RainTomorrow']
df_Dir
Valori di probabilità simili un po' per tutti i valori di wind gust direction, in particolare per quelli che procedono nella stessa direzione. Provo a discretizzarli e assegnargli valori numerici, raggruppando ogni valore nella rispettiva direzione cardinale; quindi la discretizzazione viene effettuata considerando come uniche direzioni il Nord,Sud,Est e Ovest.
df.WindGustDir.replace({'N[A-Z]*':0,'E[A-Z]*':1,'S[A-Z]*':2,'W[A-Z]*':3},regex=True,inplace=True)
df_Dir = df.groupby(['WindGustDir', 'RainTomorrow'])[['WindGustSpeed']].count()
df_Dir = df_Dir.rename({'WindGustSpeed': 'count'}, axis=1)
df_Dir['prob'] = df_Dir['count'] / df.groupby('WindGustDir').count()['RainTomorrow']
df_Dir
Anche in questo caso asumono probabilità simili, prime di vedere se influisce su WindGustSpeed discretizzo tutti gli attributi su cui è possibile farlo.
Prima di concludere le nostre analisi sul dataset, controllo la distribuzione dei valori degli attributi nello spazio. Questo può essere utile per capire se ci sono intervalli che raggrupppano un gran numero di entries, se così fosse converrebbe discretizzare gli attributi per migliorare le performance.
df.hist(bins = 50, figsize = (20,15) )
plt.show()
Distribuzioni gaussiane, fatta eccezine per RainFall che presenta un enorme picco in zero e WindGustSpeed che present vari picchi. Questi dati ci fanno capire che, molto probabilmente, se discretizziamo i dati dovremmo ottenere ottimi risultati, considerando che per alcuni valori sono presenti un numero spropositato di istanze. Prima di fare questo utilizo una pipeline nella quale insirisco al posto dei null la media, tanto per quei valori rimasti la percentuale massima di null non supera il 6%.
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
pipeline = Pipeline([
('imputer', SimpleImputer(strategy="median"))])
columns=df.columns
df=pipeline.fit_transform(df)
df=pd.DataFrame(df,columns=columns)
df.RainToday = df.RainToday.astype(int)
df.RainTomorrow = df.RainTomorrow.astype(int)
df.info()
Adesso discretizzo e verifico che dopo la discretizzazione sia tutto a posto.
df['Pressure'] = pd.cut(df.Pressure, bins = [970., 1000., 1010. ,1020.,1030.,np.inf], labels=[0,1,2,3,4])
tmp = df.groupby(['Pressure', 'RainTomorrow'])[['RainToday']].count()
tmp = tmp.rename({'RainToday':'count'}, axis=1)
tmp['prob'] = tmp['count']/df.groupby('Pressure').count()['RainTomorrow']
tmp
Ottima discretizzazione, più la pressione aumenta, più la probabilità di piovere diminuisce.
df['Humidity'] = pd.cut(df.Humidity, bins = [-0.1, 20., 40. ,60.,80.,np.inf], labels=[0,1,2,3,4])
tmp = df.groupby(['Humidity', 'RainTomorrow'])[['RainToday']].count()
tmp = tmp.rename({'RainToday':'count'}, axis=1)
tmp['prob'] = tmp['count']/df.groupby('Humidity').count()['RainTomorrow']
tmp
Qui avviene il contrario, più l'umidità è alta più è probabile che piova.
df['WindGustSpeed'] = pd.cut(df.WindGustSpeed, bins = [0., 20., 35.,40.0 ,60.,80.,np.inf], labels=[0,1,2,3,4,5])
tmp = df.groupby(['WindGustSpeed', 'RainTomorrow'])[['RainToday']].count()
tmp = tmp.rename({'RainToday':'count'}, axis=1)
tmp['prob'] = tmp['count']/df.groupby('WindGustSpeed').count()['RainTomorrow']
tmp
Più aumenta il windgustspeed più aumenta la probabilità di piovere, anche se di poco.
Adesso controllo meglio Rainfall per capire come discretizzarlo, dato che ha il 50% di valori a 0 (picco elevato).
ax=df.Rainfall.hist(bins = 50, figsize = (10,5))
ax.set_xticks([i for i in range(0,50,10) ])
plt.show(fig)
df['Rainfall'] = pd.cut(df.Rainfall, bins = [-0.1, 9.,20.,np.inf], labels=[0,1,2])
tmp = df.groupby(['Rainfall', 'RainTomorrow'])[['RainToday']].count()
tmp = tmp.rename({'RainToday':'count'}, axis=1)
tmp['prob'] = tmp['count']/df.groupby('Rainfall').count()['RainTomorrow']
tmp
Rainfall dopo averlo discretizzato, si è rivelato meglio di quanto pensassi.
df['MaxTemp'] = pd.cut(df.MaxTemp, bins = [-100., 10., 20., 34.,np.inf], labels=[0,1,2,3])
temp = df.groupby(['MaxTemp', 'RainTomorrow'])[['RainToday']].count()
temp = temp.rename({'RainToday':'count'}, axis=1)
temp['prob'] = temp['count']/df.groupby('MaxTemp').count()['RainTomorrow']
temp
Non piove se la temperatura massima è tra [34,+inf) (90%), dopo i 20° comunque, la probabilità che non piova è circa intorno all'80%.
df['MinTemp'] = pd.cut(df.MinTemp, bins = [-10., 5., 20., np.inf], labels=[0,1,2])
temp = df.groupby(['MinTemp', 'RainTomorrow'])[['RainToday']].count()
temp = temp.rename({'RainToday':'count'}, axis=1)
temp['prob'] = temp['count']/df.groupby('MinTemp').count()['RainTomorrow']
temp
Non piove se la temperatura minima è tra [-10,5], 85% di probabilità.
La temperatura non sembra molto influente sulla pioggia, ma può essere utile per gli altri attributi. Per verificare che gli attributi scelti fin'ora vanno bene, utilizzo il random forest, per capire se ognuno di essi è fondamentale.
from sklearn.ensemble import RandomForestClassifier
rnd_clf = RandomForestClassifier(n_estimators=500, random_state=42)
df_train=df.copy()
df_train=df_train.drop('RainTomorrow',axis=1)
label=df['RainTomorrow']
rnd_clf.fit(df_train,label)
nameCol=df_train.columns
fig, ax = plt.subplots(1,1)
fig.set_size_inches(9,9)
ax.plot(rnd_clf.feature_importances_,'b-')
ax.set_xticks([i for i in range(len(nameCol))])
xlabels = ax.set_xticklabels(nameCol,rotation=30, fontsize=12)
_=ax.set_title('Importanza degli attributi')
for name, score in zip(nameCol, rnd_clf.feature_importances_):
print(name, score)
Tutti gli attributi hanno egual importanza, considerando che lo scarto tra di loro è minimo. Fatta eccezione per l'attributo Cloud, che è leggermente più alto. Ma la vera eccezzione è l'attributo Humidity, il quale è significativamente più alto, e quindi migliore, degli altri. Attributo fondamentale.
Dopo aver effettuato tutte le modifiche agli attributi e dopo esserci accertati dell'effettiva importanza degli attributi rimanenti, abbiamo ottenuto la versione finale, che verrà utilizzata durante la classificazione, del nostro dataset.
df
Controlliamo anche la nuova distribuzione degli attributi.
for col in df.columns :
df[col]=df[col].values.astype('int64')
df.hist(figsize=[15,15])
Una volta ottenuti i nostri attributi di riferimento, possiamo considerare concluso il preprocessing. L'ultimo passo è quello di dividere il test set dal train set, mantenendo in percentuale la stessa distribuzione di valori, considerando l'attributo sul quale effettueremo la classificazione, ovvero RainTomorrow. Con scikitleran divido il dataset in train e test e controllo se la divisione rispecchia i parametri citati qui sopra.
from sklearn.model_selection import train_test_split
training_set, test_set = train_test_split(df, test_size=0.2, random_state=42)
training_set["RainTomorrow"].hist()
test_set["RainTomorrow"].hist()
Split perfetto! Adesso creo le label
training_label=training_set['RainTomorrow'].copy()
training_set=training_set.drop("RainTomorrow",axis=1)
test_label=test_set['RainTomorrow'].copy()
test_set=test_set.drop('RainTomorrow',axis=1)
Ora che abbiamo un dataset utile per i nostri scopi, non ci resta che scoprire quale modello di classificazione riesce a predire correttamente, dato un record, se il giorno dopo pioverà oppure no.
Questa è una sezione extra, non era richiesto e non era necessiario effettuare del clustering. Però, per completezza abbiamo deciso di utilizzare anche questo stumento.
Si è provato a fare una analisi per località, l'obbiettivo era quello di trovare zone geografiche che avessero una piovosità simile, in realtà non si è ottenuto il risultato richiesto, probabilmente perchè i dati relativi a data e località sono abbastanza incompleti, come già mostrato nella fase di preprocessing.
Dunque il clustering viene fatto sui dati ottenuti nella fase del preprocessing.
# clus lo utilizzo per il clustering
clus = training_set.copy()
Il primo algoritmo utilizzato per il clustering è il DBSCAN, si è pensato che, parlando di dati spazio-temporali, esso potesse modellare abbastanza bene i punti che dovrebbero avere le stesse densità.
from sklearn.cluster import DBSCAN # dbscan
clusters = np.array([])
for e in range(7,11):
ep = e/10
dbscan = DBSCAN(eps = ep, min_samples=150)
dbscan.fit(clus)
labels = pd.Series(dbscan.labels_)
labels = pd.Series(labels.unique())
labels = np.ones_like(labels)
clusters = np.append(clusters, labels.sum())
clusters
Il numero di classi è ciò che ci aspettiamo. Avendo un dataset molto denso, è normale aspettarsi che per un eps < 1 abbiamo un numero di cluster (65) costante indipendente dal variare di eps, mentre per eps >= 1 succede la stessa cosa, ma il numero di cluster è nettamente più basso (2) proprio perchè i punti sono molto densi.
Vediamo allora cosa succede facendo variare min_samples.
clusters = np.array([])
for e in range(150, 300, 50):
dbscan = DBSCAN(eps = 0.8, min_samples=e)
dbscan.fit(clus)
labels = pd.Series(dbscan.labels_)
labels = pd.Series(labels.unique())
labels = np.ones_like(labels)
clusters = np.append(clusters, labels.sum())
clusters
plt.figure(figsize=(8, 3.5))
plt.plot(range(3), clusters, "bo-")
plt.xlabel("$Min Samples$", fontsize=14)
plt.ylabel("Cluster Numbers", fontsize=14)
y_min = clusters.min()- 5
y_max = clusters.max()+ 5
plt.axis([-0.5, 2.5, y_min, y_max ])
plt.show()
clusters[1]
Scelgo come eps valida 0.8, che mi permette di avere un buon numero di cluster (39), quando min_samples = 200, si è scelto di avere questo numero di cluster in base alla correttezza di Kmeans che viene eseguito dopo, in modo da avere un numero di cluster che ipoteticamente migliorano l'errore il più possibile.
dbscan = DBSCAN(eps = 0.8, min_samples=200)
clusters = dbscan.fit_predict(clus)
cl = pd.Series(clusters)
labels = cl.unique()
labels
cl.hist(bins=40, figsize=[10,10])
Ci sono 39 cluster. I punti mostrano uno sbilanciamento netto, questo è sicuramente dato dalla forte densità dei punti (come ci si aspetta). Gli altri adattamenti provati mostrano comunque un fortissimo sbilanciamento, ma con un numero di cluster terribilemtne più alto abbassando min_samples, oppure terribilmente più basso portando eps>=1, dunque questo dovrebbe essere un adattamento, tutto sommato, buono.
Passiamo a KMeans e andiamo a valutare i paramentri in input in base alla discesa dell'SSE medio.
from sklearn.cluster import KMeans
inertia = np.array([])
for k in range(36,42):
kmeans =KMeans(k).fit(clus)
inertia = np.append(inertia, kmeans.inertia_)
inertia
plt.figure(figsize=(8, 3.5))
plt.plot(range(36, 42), inertia, "bo-")
plt.xlabel("$k$", fontsize=14)
plt.ylabel("Inertia", fontsize=14)
y_max = inertia.max() + 10000
y_min = inertia.min() - 10000
plt.axis([35.5, 41.5, y_min, y_max])
plt.show()
L'algoritmo con valori di K più grandi fa fatica a terminare, quindi mi sono mantenuto su un valore di cluster che vanno da 45 a 51, dopo diverse run ho notato come a k=49 l'SSE medio mostra sempre un buon abbassamento.
Anche se comunque l'errore è alto, non possiamo fare di molto meglio e quindi manteniamo k=49.
kmeans = KMeans(40).fit(clus)
centroidi = kmeans.cluster_centers_
def plot_centroids(centroids, weights=None, circle_color='w', cross_color='k'):
if weights is not None:
centroids = centroids[weights > weights.max() / 10]
plt.scatter(centroids[:, 0], centroids[:, 1],
marker='o', s=30, linewidths=8,
color=circle_color, zorder=10, alpha=0.9)
plt.scatter(centroids[:, 0], centroids[:, 1],
marker='x', s=50, linewidths=50,
color=cross_color, zorder=11, alpha=1)
plot_centroids(centroidi)
kmeans = kmeans.predict(clus)
kmeans
cl = pd.Series(kmeans)
cl.unique()
cl.hist(bins=40, figsize=[10,10])
Sicuramente i cluster sono bilanciati utilizzando kmeans, a volte qualche cluster è vuoto, segno che i dati sono molto densi.
Non sembrano esserci grandi risultati dalla creazione dei cluster (sia con DBSCAN, che con KMeans), tuttavia proveremo questi adattamenti con alcuni algoritmi di classificazione per valutare se effettivamnete non abbiamo miglioramenti.
from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
def compute_performance(estimators, X, y, scoring_funs, cross_val=True):
'''
Compute every score function for each estimator
Parameters:
----------
estimators : a list of predictors
X : data used for computed the corss_validation_test
y : target values
scoring_fun : list of performances measures
Returns:
--------
A dictionary containing the score of each estimator
'''
score_dict = {}
for e in estimators:
score_dict[e.__class__.__name__] = {}
if not cross_val:
y_pred = e.predict(X)
for f in scoring_funs:
if cross_val:
score = cross_val_score(e, X ,y, cv=10 ,scoring=f).mean()
name = f
else:
_f, name = f
score = _f(y_pred, y)
score_dict[e.__class__.__name__][name] = score
return score_dict
Il primo classificatore che andiamo ad addestrare è l'albero decisionale utilizzando la metrica del GinIndex.
from sklearn.tree import DecisionTreeClassifier
DTree_clf = DecisionTreeClassifier().fit(training_set, training_label)
nameCol=training_set.columns
fig, ax = plt.subplots(1,1)
fig.set_size_inches(9,9)
ax.plot(DTree_clf.feature_importances_,'b-')
ax.set_xticks([i for i in range(len(nameCol))])
xlabels = ax.set_xticklabels(nameCol,rotation=30, fontsize=12)
_=ax.set_title('Importanza degli attributi')
Confermata l'importanza degli attributi trovata in random forest.
impurità = pd.Series(DTree_clf.tree_.impurity)
print("Frazione di nodi puri: %.2f %%\n" %((1 - impurità[impurità > 0.].count()/impurità.count())*100))
Poco più del 26 % dei nodi dell'albero sono puri.
compute_performance([DTree_clf], training_set, training_label, scoring_funs=['accuracy', 'recall', 'precision', 'f1'])
test_score = DTree_clf.score(test_set,test_label)
training_score = DTree_clf.score(training_set,training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
Misurando lo score ottenuto sul training set e sul test set, ci rendiamo conto di come il modello si adatti bene ad entrambi, per questo motivo supponiamo che non ci sia overfitting.
Ora proveremo ad applicare il modello clusterizzando i dati con l'algoritmo KMeans.
pipeline = Pipeline([
("kmeans", KMeans(n_clusters=49)),
("Decision_tree", DecisionTreeClassifier()),
])
pipeline.fit(training_set, training_label)
test_score = pipeline.score(test_set,test_label)
training_score = pipeline.score(training_set,training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
DTree_clf = pipeline.named_steps.Decision_tree
compute_performance([DTree_clf], training_set, training_label, scoring_funs=['accuracy', 'recall', 'precision', 'f1'])
Non ci sono varianzioni di performance rispetto al modello precedente, dunque la clusterizzazione non ha effetto.
Proviamo adesso ad untilizzare l'entropia invece del GiniIndex.
DTree_clf = DecisionTreeClassifier(criterion='entropy').fit(training_set, training_label)
test_score = DTree_clf.score(test_set,test_label)
training_score = DTree_clf.score(training_set,training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
compute_performance([DTree_clf], training_set, training_label, scoring_funs=['accuracy', 'recall', 'precision', 'f1'])
impurità = pd.Series(DTree_clf.tree_.impurity)
print("Frazione di nodi puri: %.2f %%\n" %((1 - impurità[impurità > 0.].count()/impurità.count())*100))
Anche il numero di nodi puri è similare.
Il classificatore che utilizzaimo adesso implementa modelli lineari regolarizzati con apprendimento discendente del gradiente stocastico. Il modello è della familgia dei classificatori Support Vector Machines.
from sklearn.linear_model import SGDClassifier
sgd_clf = SGDClassifier(max_iter=1000, tol=-np.infty)
sgd_clf.fit(training_set, training_label)
compute_performance([sgd_clf], training_set, training_label, ['accuracy', 'precision', 'recall', 'f1'])
Miglioriamo di molto sulla precision e scendiamo di poco sulla recall.
coef = pd.DataFrame({'Frequenza Coefficenti SGDC':sgd_clf.coef_[0]})
coef.hist()
test_score = sgd_clf.score(test_set,test_label)
training_score = sgd_clf.score(training_set,training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
Anche qua non sembra esserci overfitting, ma c'è un peggioramento di score rispetto all'albero decisionale.
Proviamo ad addestrare un classificatore Bayesiano.
from sklearn.naive_bayes import GaussianNB
nb_clf = GaussianNB()
nb_clf = nb_clf.fit(training_set, training_label)
compute_performance([nb_clf], training_set, training_label, ['accuracy', 'precision', 'recall', 'f1'])
Risultati non molto buoni dal punto di vista della precision (scendiamo sia rispetto al DTree che al SGD), ma milgioriamo di qualche punto rispetto alla recall.
test_score = nb_clf.score(test_set,test_label)
training_score = nb_clf.score(training_set,training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
Lo score somiglia di più a quello ottenuto con SGD che a quello del Dtree, la differenza tra training set e test set questa volta è minima.
Prima di Addestrare il classificatore basato sui KNN, si è applicata una fase pre addestramento volta a cercare il miglior adattamento in termini di numero di knn.
from sklearn.neighbors import KNeighborsClassifier
scores = []
for k in range(10, 100, 20):
knn_clf = KNeighborsClassifier(n_neighbors=k)
knn_clf.fit(training_set, training_label)
scores.append(cross_val_score(knn_clf, training_set, training_label, cv=10 ,scoring='recall').mean())
scores
Si è voluto misurare la Recall, essendo esso il parametro che più si dovrebbe migliorare, tuttavia cambiare il numero di KNN in input non dà grandi vantaggi in questo senso. si è scelto KNN=60.
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier(n_neighbors=60)
knn_clf.fit(training_set, training_label)
compute_performance([knn_clf], training_set, training_label, ['accuracy', 'precision', 'recall', 'f1'])
test_score = knn_clf.score(test_set,test_label)
training_score = knn_clf.score(training_set,training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
Lo score è abbastanza alto (più dei due precedenti, ma sempre meno del DTree) con una differenza tra training set e test set, minima.
Proviamo a Clusterizzare con DBSCAN prima di addestrare il modello.
dbscan = DBSCAN(eps= 0.8, min_samples=200)
dbscan = dbscan.fit(training_set)
knn = KNeighborsClassifier(n_neighbors=60)
knn.fit(dbscan.components_, dbscan.labels_[dbscan.core_sample_indices_])
compute_performance([knn], training_set, training_label, ['accuracy', 'precision', 'recall', 'f1'])
test_score = knn.score(test_set,test_label)
training_score = knn.score(training_set,training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
In temrini di Accuracy, Precision, Recall e F-measure non ci sono cambiamenti, oltretutto lo score degrada moltissimo, sicuramente non un buon risultato.
Proviamo ad addestrare il classico classificatore Support Vector.
from sklearn.svm import LinearSVC
svc_clf = LinearSVC(max_iter= 500)
svc_clf.fit(training_set, training_label)
compute_performance([svc_clf], training_set, training_label, ['accuracy', 'precision', 'recall', 'f1'])
test_score = svc_clf.score(test_set,test_label)
training_score = svc_clf.score(training_set,training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
Risultati comunque abbastanza buoni, meglio delgi altri classificatori in generale.
Cerchiamo di migliorare i risultati andando ad utilizzare gli ensemble. Nello specifico verranno addestrati: RandomForest, Voting e Bagging.
from sklearn.ensemble import RandomForestClassifier
rnd_clf = RandomForestClassifier(n_estimators=500)
rnd_clf.fit(training_set,training_label)
compute_performance([rnd_clf], training_set, training_label, scoring_funs=['accuracy', 'recall', 'precision'])
Le statistiche del RandomForest sono esattamente uguali a quelle del semplice DTree.
Per il Voting classifier, si sono utilizzati i tre classificatori visti prima: SVC, NaiveBayes e KNN, e in più si è aggiunto anche il Random Forest.
from sklearn.ensemble import VotingClassifier
nb_clf = GaussianNB()
svc_clf = LinearSVC(max_iter= 500)
rnd_clf = RandomForestClassifier(n_estimators=500)
knn_clf = KNeighborsClassifier(n_neighbors=60)
voting_clf = VotingClassifier(
estimators=[('nb', nb_clf), ('svc', svc_clf), ('rf', rnd_clf), ('knn', knn_clf)],
voting='hard')
voting_clf.fit(training_set,training_label)
compute_performance([voting_clf], training_set, training_label, scoring_funs=['accuracy', 'recall', 'precision'])
Il migliore in termini di precision, ma perdiamo troppo sulla recall.
Infine il BaggingClassifier viene utilizzato con il Dtree, si è provato anche con altri modelli, ma non si sono ottenuti miglioramenti.
from sklearn.ensemble import BaggingClassifier
bag_clf = BaggingClassifier(
DecisionTreeClassifier(), n_estimators=500,
max_samples=100, bootstrap=True)
bag_clf.fit(training_set, training_label)
compute_performance([bag_clf], training_set, training_label, scoring_funs=['accuracy', 'recall', 'precision'])
Il Bagging è quello che ha dato il miglior risultato, ripettando le caratteristiche dei dati.
Il problema principale dei classificatori sono sulle recall, molto basso come valore. Arrivati a aquesto punto si suppone che non si possa fare molto meglio di così.
Adesso proviamo a calcolare matrice di confusione e curva ROC di quest'ultimo modello.
from sklearn.metrics import confusion_matrix
from sklearn.utils.multiclass import unique_labels
def plot_confusion_matrix(y_true, y_pred, classes,
classes_value,
normalize=False,
title=None,
cmap=plt.cm.Blues):
"""
This function prints and plots the confusion matrix.
Normalization can be applied by setting `normalize=True`.
"""
if not title:
if normalize:
title = 'Normalized confusion matrix'
else:
title = 'Confusion matrix, without normalization'
# Compute confusion matrix
cm = confusion_matrix(y_true, y_pred)
# Only use the labels that appear in the data
classes = classes[unique_labels(y_true, y_pred)]
if normalize:
cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
print("Normalized confusion matrix")
else:
print('Confusion matrix, without normalization')
print(cm)
fig, ax = plt.subplots()
im = ax.imshow(cm, interpolation='nearest', cmap=cmap)
ax.figure.colorbar(im, ax=ax)
# We want to show all ticks...
ax.set(xticks=np.arange(cm.shape[1]),
yticks=np.arange(cm.shape[0]),
# ... and label them with the respective list entries
xticklabels=classes_value, yticklabels=classes_value,
title=title,
ylabel='True label',
xlabel='Predicted label')
# Rotate the tick labels and set their alignment.
plt.setp(ax.get_xticklabels(), rotation=45, ha="right",
rotation_mode="anchor")
# Loop over data dimensions and create text annotations.
fmt = '.2f' if normalize else 'd'
thresh = cm.max() / 2.
for i in range(cm.shape[0]):
for j in range(cm.shape[1]):
ax.text(j, i, format(cm[i, j], fmt),
ha="center", va="center",
color="white" if cm[i, j] > thresh else "black")
fig.tight_layout()
return ax
classes_value = ["Rain Tomorrow: NO", "RainTomorrow: YES"]
test_pred = bag_clf.predict(test_set)
plot_confusion_matrix(test_label, test_pred, classes=df.RainTomorrow, classes_value=classes_value)
plot_confusion_matrix(test_label, test_pred, classes=df.RainTomorrow, classes_value=classes_value, normalize=True)
from sklearn.metrics import roc_curve
def plot_roc_curve(fpr, tpr, label=None):
plt.plot(fpr, tpr, linewidth=2, label=label)
plt.plot([0, 1], [0, 1], 'k--')
plt.axis([0, 1, 0, 1])
plt.xlabel('False Positive Rate', fontsize=16)
plt.ylabel('True Positive Rate', fontsize=16)
from sklearn.model_selection import cross_val_predict
train_scores = cross_val_predict(bag_clf, training_set, training_label, cv=3)
fpr, tpr, thresholds = roc_curve(training_label, train_scores)
plt.figure(figsize=(8, 6))
plot_roc_curve(fpr, tpr)
plt.show()
Le prestazioni non sono molto soddisfacenti, per non parlare della precision e, soprattutto, la recall che sono molto basse.Quindi abbiamo decisio di effettuare delle operazioni di post-processing. Per migliorare questi fattori abbiamo provato a tornare indetro al dataset originale, e a compiere delle scelti differenti er quanto riguarda gli attributi. Nonostante le nostre innumerevoli scelte non si sono riscontrati dei miglioramenti, anzi in molti casi peggioravano le prestazioni. L'ultima nostra possibilità per migliorare le prestuazioni, riguradava il bilanciamento del dataset. Perchè sii è notato che nel dataset origninale c'è un forte sbilanciamento delle classi (80% di No e 20% di Si), probabilmente a causa di questo la classe 'Yes' viene vista come classe rara anche se non dovrebbe. Dunque proviamo a creare un dataset che sia più bilanciato (40 % di si) e vediamo se riusciamo ad ottenere buone performance
Abbiamo creato una funzione, che a partire dal dataset originale, crea un dataset bilanciato secondo le nostre esigenze.
from sklearn.utils import shuffle
def creaDatasetBilanciato(df, class_coloumn, num_elem=5000,bilancio=0.35):
df_yes = df[df[class_coloumn] == 1]
df_no = df[df[class_coloumn] == 0]
sample_yes = int(num_elem*bilancio)
sample_no = int(num_elem*(1-bilancio))
df_yes = df_yes.sample(n=sample_yes)
df_no = df_no.sample(n=sample_no)
df_new = pd.concat([df_yes, df_no])
df_new = shuffle(df_new, random_state=1)
return df_new
df_ = creaDatasetBilanciato(df, num_elem=90000, class_coloumn='RainTomorrow')
df_.hist(figsize=[15,15])
La distribuzione degli attributi è uguale a quella del dataset originale, quindi la funzione è fatta bene.
_training_set, _test_set = train_test_split(df_, test_size=0.2, random_state=42)
_training_set["RainTomorrow"].hist()
_test_set["RainTomorrow"].hist()
_training_label=_training_set['RainTomorrow'].copy()
_training_set=_training_set.drop("RainTomorrow",axis=1)
_test_label=_test_set['RainTomorrow'].copy()
_test_set=_test_set.drop('RainTomorrow',axis=1)
DTree_clf = DecisionTreeClassifier().fit(_training_set, _training_label)
compute_performance([DTree_clf], _training_set, _training_label, scoring_funs=['accuracy', 'recall', 'precision', 'f1'])
impurità = pd.Series(DTree_clf.tree_.impurity)
print("Frazione di nodi puri: %.2f %%\n" %((1 - impurità[impurità > 0.].count()/impurità.count())*100))
test_score = DTree_clf.score(_test_set,_test_label)
training_score = DTree_clf.score(_training_set,_training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
Si abbassa di molto l'accuracy (circa del 12-13%), però c'è un netto milgioramento per quanto riguarda la recall che passa dal 40% al 56% perdendo solo il 3% sulla precision. Si alza anche la percentuale di nodi puri.
Tuttavia lo score si abbassa (così coma l'accuracy), questo può essere sintomo di un aumento di overfitting, anche se in realtà abbiamo notato che aumentando il numero di record (num_ele) nella creazione del dataset bilanciato si migliora questo aspetto. Tuttavia non possiamo alzare il parametro perchè non abbiamo a disposizione troppi record con etichetta di classe 'YES'.
sgd_clf = SGDClassifier(max_iter=1000, tol=-np.infty)
sgd_clf.fit(_training_set, _training_label)
compute_performance([sgd_clf], _training_set, _training_label, ['accuracy', 'precision', 'recall', 'f1'])
test_score = sgd_clf.score(_test_set,_test_label)
training_score = sgd_clf.score(_training_set,_training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
nb_clf = GaussianNB()
nb_clf = nb_clf.fit(_training_set, _training_label)
compute_performance([nb_clf], _training_set, _training_label, ['accuracy', 'precision', 'recall', 'f1'])
test_score = nb_clf.score(_test_set,_test_label)
training_score = nb_clf.score(_training_set,_training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
knn_clf = KNeighborsClassifier(n_neighbors=60)
knn_clf.fit(_training_set, _training_label)
compute_performance([knn_clf], _training_set, _training_label, ['accuracy', 'precision', 'recall', 'f1'])
test_score = knn_clf.score(_test_set,_test_label)
training_score = knn_clf.score(_training_set,_training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
svc_clf = LinearSVC(max_iter= 500)
svc_clf.fit(_training_set, _training_label)
compute_performance([svc_clf], _training_set, _training_label, ['accuracy', 'precision', 'recall', 'f1'])
test_score = svc_clf.score(_test_set,_test_label)
training_score = svc_clf.score(_training_set,_training_label)
print("Training Score : %.3f ," %training_score, " Test Score : %.3f " %test_score)
Proviamo anche gli ensamble learnig su questo nuovo training set.
rnd_clf = RandomForestClassifier(n_estimators=500)
rnd_clf.fit(_training_set,_training_label)
compute_performance([rnd_clf], _training_set, _training_label, scoring_funs=['accuracy', 'recall', 'precision'])
nb_clf = GaussianNB()
svc_clf = LinearSVC(max_iter= 500)
rnd_clf = RandomForestClassifier(n_estimators=500)
knn_clf = KNeighborsClassifier(n_neighbors=60)
voting_clf = VotingClassifier(
estimators=[('nb', nb_clf), ('svc', svc_clf), ('rf', rnd_clf), ('knn', knn_clf)],
voting='hard')
voting_clf.fit(_training_set,_training_label)
compute_performance([voting_clf], _training_set, _training_label, scoring_funs=['accuracy', 'recall', 'precision'])
bag_clf = BaggingClassifier(
DecisionTreeClassifier(), n_estimators=500,
max_samples=100, bootstrap=True)
bag_clf.fit(_training_set, _training_label)
compute_performance([bag_clf], _training_set, _training_label, scoring_funs=['accuracy', 'recall', 'precision'])
classes_value = ["Rain Tomorrow: NO", "RainTomorrow: YES"]
test_pred = bag_clf.predict(_test_set)
plot_confusion_matrix(_test_label, test_pred, classes=df.RainTomorrow, classes_value=classes_value)
plot_confusion_matrix(_test_label, test_pred, classes=df.RainTomorrow, classes_value=classes_value, normalize=True)
Come prevedebile sale di poco la percentuale di FN e scende abbastanza quella dei FP, allo stesso modo scende la percentuale di TP e sale quella di TN
train_scores = cross_val_predict(bag_clf, _training_set, _training_label, cv=3)
fpr, tpr, thresholds = roc_curve(_training_label, train_scores)
plt.figure(figsize=(8, 6))
plot_roc_curve(fpr, tpr)
plt.show()
Il punto di picco si sposta da 0.4 a 0.6, mentre l'andamento resta simile a prima.
Anche questa sezione si potrebbe considerare extra dato che non era obbligatoria. Però le reti neurali in genere forniscono degli ottimi risultati come classificatori, quindi abbiamo voluto testare diverse reti per controllare se le performance aumentano.
La rete neurale presente qui sotto è formata da 3 strati: lo strato di input, intermedio e di output. Ho provato anche delle reti di diversa tipologia, deep leraning ad esempio, però questa è quella che presenta le performance migliori. I primi due strati sono formati da quattro neuroni e la loro funzione di attivazione è la 'relu', è quella che presentava i risultati migliori. Mentre, l'ultimo strato, quello che fornisce i risultati della rete, è formato da due neuroni che rappresentano rispettivamente la classe 'no' e la classe 'yes'. Inoltre grazie alla funzione di attivazione 'softmax' oltre alla classe, restituisce in output la probabilità associatà alla classe. Quindi è come se fornisse una misura che indica la sicurezza che ha la rete nell'assegnare la classe a un record.
from keras.utils import to_categorical
test_label2=test_label.copy()
#one-hot encode labels
training_label=to_categorical(training_label)
test_label=to_categorical(test_label)
import keras
modello=keras.Sequential()
modello.add(keras.layers.Dense(units=4,activation='relu',input_shape=(9,)))
modello.add(keras.layers.Dense(units=4,activation='relu'))
modello.add(keras.layers.Dense(units=2,activation='softmax'))
epoche=100
learning_rate=0.01
batch=1024
modello.compile(optimizer=keras.optimizers.SGD(learning_rate),
loss='categorical_crossentropy',
metrics=['accuracy']
)
history=modello.fit(training_set,training_label, validation_data=(test_set, test_label),epochs=epoche,batch_size=batch)
plt.figure(num=None, figsize=(8, 6), dpi=80)
plt.title('Funzione di Costo')
plt.plot(history.history['loss'], label='train')
plt.plot(history.history['val_loss'], label='test')
plt.xlabel('Epoca')
plt.legend()
plt.show()
plt.figure(num=None, figsize=(8, 6), dpi=80)
plt.title('Accuracy')
plt.plot(history.history['acc'], label='train')
plt.plot(history.history['val_acc'], label='test')
plt.xlabel('Epoca')
plt.legend()
plt.show()
from sklearn.metrics import accuracy_score
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import f1_score
y_classes = modello.predict_classes(test_set, verbose=0)
accuracy = accuracy_score(test_label2, y_classes)
print('Accuracy: %f' % accuracy)
precision = precision_score(test_label2, y_classes)
print('Precision: %f' % precision)
recall = recall_score(test_label2, y_classes)
print('Recall: %f' % recall)
f1 = f1_score(test_label2, y_classes)
print('F1 score: %f' % f1)
Risultati in linea con quelli ottenuti con gli altri classificatori, solamente la recall migliora leggermente.
from sklearn.metrics import precision_recall_curve
precisions, recalls, thresholds = precision_recall_curve(test_label2, y_classes)
plt.figure(figsize=(8, 6))
plt.plot(recalls, precisions, "b-", linewidth=2)
plt.xlabel("Recall", fontsize=16)
plt.ylabel("Precision", fontsize=16)
plt.axis([0, 1, 0, 1])
plt.show()
Notiamo da questa curva, come sia difficile ottenere una buona recall e una buona precision allo stesso tempo, è come se fossere inversamente proporzionali.
from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(test_label2, y_classes)
plt.figure(figsize=(8, 6))
plt.plot(fpr, tpr, linewidth=2, label=label)
plt.plot([0, 1], [0, 1], 'k--')
plt.axis([0, 1, 0, 1])
plt.xlabel('False Positive Rate', fontsize=16)
plt.ylabel('True Positive Rate', fontsize=16)
plt.show()
Esattamente la stessa curva ROC presente nel miglior modello della classificazione.
from sklearn.metrics import confusion_matrix
import itertools
def plot_confusion_matrix(cm, classes, normalize=False, title='Confusion matrix', cmap=plt.cm.Blues):
plt.imshow(cm, interpolation='nearest', cmap=cmap)
plt.title(title)
plt.colorbar()
tick_marks = np.arange(len(classes))
plt.xticks(tick_marks, classes, rotation=45)
plt.yticks(tick_marks, classes)
fmt ='d'
thresh = cm.max() / 2.
for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
plt.text(j, i, format(cm[i, j], fmt),
horizontalalignment="center",
color="white" if cm[i, j] > thresh else "black")
plt.tight_layout()
plt.ylabel('True label')
plt.xlabel('Predicted label')
plt.show()
confusion = confusion_matrix(test_label2, y_classes)
plot_confusion_matrix(confusion, classes=['Rain Tomorrow: NO', 'Rain Tomorrow: YES'], title='Confusion matrix')
Anche la matrice di confusione, ovviamente presenta dei risultati simili a quelli ottenuti con la classificazione.
Il vantaggio della rete neurale è che oltre alla classe fornisce la probabilità di pioggia, quindi un'utente può decidere, grazie al valore della probabilità, se fidarsi o meno di una determinata classificazione. Ecco di seguito un'esempio di output fornito dalla rete:
output=["No Rain Tomorrow","Rain Tomorrow"]
classe=modello.predict_classes(test_set[0:5], verbose=0)
prob=modello.predict(test_set[0:5],verbose=0)
for p,c in zip(prob,classe):
print('Classificazione-> Classe: '+output[c]+', Probabilità: %.1f%s' % (p[c]*100,'%'))
L'utilizzo della rete neurale non ha portato a delle migliorie rilevanti, per questo abbiamo eseguito le stesse operazioni di post-processing effettuate con la classificazione. Però, oltre all'utilizzo del dataset bilanciato, ho effettuato un'operazione sulle predizioni effettuate sul dataset orignale. Siccome la rete fornisce anche la probabilità e abbiamo problemi con la Recall, quindi con i falsi negativi, ho deciso di cambiare la classe 'no' in 'yes' qunado la rete è poco sicura. Ovvero quando la probabilità di 'no' è compresa tra il 50% e 60% cambia la classe in 'yes', perchè vuol dire che non è molto sicura su come classificare. Una volta fatta quest'operazione ho misurato le performance.
ymod_classes = modello.predict_classes(test_set, verbose=0)
ymod_probs = modello.predict(test_set, verbose=0)
index=0
for classe,prob in zip(ymod_classes,ymod_probs) :
if classe==0 and (0.50<prob[0]<0.59) :
ymod_classes[index]=1
ymod_probs[index,1]=1-ymod_probs[index,1]
index+=1
accuracy = accuracy_score(test_label2, ymod_classes)
print('Accuracy: %f' % accuracy)
precision = precision_score(test_label2, ymod_classes)
print('Precision: %f' % precision)
recall = recall_score(test_label2, ymod_classes)
print('Recall: %f' % recall)
f1 = f1_score(test_label2, ymod_classes)
print('F1 score: %f' % f1)
confusion = confusion_matrix(test_label2, ymod_classes)
plot_confusion_matrix(confusion, classes=['No Rain Tomorrow', 'Rain Tomorrow'], title='Confusion matrix')
Come previsto la recall aumenta quasi del 10% a discapito della precisione che perde circa il 5-6% in genere, mentre l'accuraci diminuisce di poco. Quindi questo metodo è molto efficace se si cerca una recall buona, quindi se all'utente interessa che vengano classificate bene le istanze della classe 'yes'.
Il resto del post-processing prosegue con il dataset bilanciato visto già nella classificazione.
from sklearn.utils import shuffle
def creaDatasetBilanciato(df, class_coloumn, num_elem=50000,bilancio=0.40):
df_yes = df[df[class_coloumn] == 1]
df_no = df[df[class_coloumn] == 0]
sample_yes = int(num_elem*bilancio)
sample_no = int(num_elem*(1-bilancio))
if df_yes.shape[0]<sample_yes:
return None
df_yes = df_yes.sample(n=sample_yes)
df_no = df_no.sample(n=sample_no)
df_new = pd.concat([df_yes, df_no])
df_new = shuffle(df_new, random_state=1)
return df_new
new_df=creaDatasetBilanciato(df,'RainTomorrow',num_elem=90000,bilancio=0.35)
if not isinstance(new_df,pd.DataFrame) :
print("Non proseguire, cambia i parametri!")
training_set, test_set = train_test_split(new_df, test_size=0.2, random_state=42)
training_set["RainTomorrow"].hist()
test_set["RainTomorrow"].hist()
training_label=training_set['RainTomorrow'].copy()
training_set=training_set.drop("RainTomorrow",axis=1)
test_label=test_set['RainTomorrow'].copy()
test_set=test_set.drop('RainTomorrow',axis=1)
test_label2=test_label.copy()
#one-hot encode labels
training_label=to_categorical(training_label)
test_label=to_categorical(test_label)
La rete è uguale, ovviamente, a quella usata in precedenza. L'unica differenza è che ora utilizzo la regolarizzazione di tipo l2 per evitare l'ioverfitting. Dato che adesso ho meno dati e quindi la rete "fatica" un po' di più ad apprendere.
modello=keras.Sequential()
parametroRegolarizzazione=0.0105
modello.add(keras.layers.Dense(units=4,kernel_regularizer=keras.regularizers.l2(parametroRegolarizzazione),activation='relu',input_shape=(9,)))
modello.add(keras.layers.Dense(units=4,kernel_regularizer=keras.regularizers.l2(parametroRegolarizzazione),activation='relu'))
modello.add(keras.layers.Dense(units=2,activation='softmax'))
epoche=175
learning_rate=0.01
batch=1024
modello.compile(optimizer=keras.optimizers.SGD(learning_rate),
loss='categorical_crossentropy',
metrics=['accuracy']
)
history=modello.fit(training_set,training_label, validation_data=(test_set, test_label),epochs=epoche,batch_size=batch)
plt.figure(num=None, figsize=(8, 6), dpi=80)
plt.title('Funzione di Costo')
plt.plot(history.history['loss'], label='train')
plt.plot(history.history['val_loss'], label='test')
plt.xlabel('Epoca')
plt.legend()
plt.show()
plt.figure(num=None, figsize=(8, 6), dpi=80)
plt.title('Accuracy')
plt.plot(history.history['acc'], label='train')
plt.plot(history.history['val_acc'], label='test')
plt.xlabel('Epoca')
plt.legend()
plt.show()
y_classes = modello.predict_classes(test_set, verbose=0)
accuracy = accuracy_score(test_label2, y_classes)
print('Accuracy: %f' % accuracy)
precision = precision_score(test_label2, y_classes)
print('Precision: %f' % precision)
recall = recall_score(test_label2, y_classes)
print('Recall: %f' % recall)
f1 = f1_score(test_label2, y_classes)
print('F1 score: %f' % f1)
precisions, recalls, thresholds = precision_recall_curve(test_label2, y_classes)
plt.figure(figsize=(8, 6))
plt.plot(recalls, precisions, "b-", linewidth=2)
plt.xlabel("Recall", fontsize=16)
plt.ylabel("Precision", fontsize=16)
plt.axis([0, 1, 0, 1])
plt.show()
fpr, tpr, thresholds = roc_curve(test_label2, y_classes)
plt.figure(figsize=(8, 6))
plt.plot(fpr, tpr, linewidth=2, label=label)
plt.plot([0, 1], [0, 1], 'k--')
plt.axis([0, 1, 0, 1])
plt.xlabel('False Positive Rate', fontsize=16)
plt.ylabel('True Positive Rate', fontsize=16)
plt.show()
confusion = confusion_matrix(test_label2, y_classes)
plot_confusion_matrix(confusion, classes=['No Rain Tomorrow', 'Rain Tomorrow'], title='Confusion matrix')
output=["No Rain Tomorrow","Rain Tomorrow"]
classe=modello.predict_classes(test_set[0:5], verbose=0)
prob=modello.predict(test_set[0:5],verbose=0)
for p,c in zip(prob,classe):
print('Classificazione-> Classe: '+output[c]+', Probabilità: %.1f%s' % (p[c]*100,'%'))
ymod_classes = modello.predict_classes(test_set, verbose=0)
ymod_probs = modello.predict(test_set, verbose=0)
index=0
for classe,prob in zip(ymod_classes,ymod_probs) :
if classe==0 and (0.50<prob[0]<0.59) :
ymod_classes[index]=1
ymod_probs[index,1]=1-ymod_probs[index,1]
index+=1
accuracy = accuracy_score(test_label2, ymod_classes)
print('Accuracy: %f' % accuracy)
precision = precision_score(test_label2, ymod_classes)
print('Precision: %f' % precision)
recall = recall_score(test_label2, ymod_classes)
print('Recall: %f' % recall)
f1 = f1_score(test_label2, ymod_classes)
print('F1 score: %f' % f1)
confusion = confusion_matrix(test_label2, ymod_classes)
plot_confusion_matrix(confusion, classes=['No Rain Tomorrow', 'Rain Tomorrow'], title='Confusion matrix')
Nell’analizzare il dataset ci si è soffermati molto nella fase di preprocessing. Oltre ad analizzare il comportamento e la correlazione di ogni attributo rispetto all’attributo di classe, vengono trasformati gli attributi categorici e ordinali in float applicando varie tecniche offerte dalle librerie di python. Alcuni attributi vengono eliminati perché hanno troppi valori null e di conseguenza non incidono nella classificazione, mentre altri sono fortemente correlati, ad esempio temp9am e minTemp mostrano quasi gli stessi valori e quindi temp9am è stata eliminata. Quest’ultimi sono dati dal fatto che probabilmente i risultati di una rilevazione temporale vengono memorizzate in entrambe le colonne (la temperatura minima è anche quella rilevata alle 9 di mattina). Dopo aver finito il preprocessing si è iniziato con la classificazione, tuttavia si è provato a applicare diverse scelte sugli attributi per provare a migliorare i parametri sulla classificazione, di queste, l’unico cambiamento che ha davvero migliorato lo score è stato la creazione di un dataset bilanciato. Quest’ultimo cambiamento si è applicato perché nel dataset originale c’è un altissimo numero di record con etichetta ‘no’ e pochi con etichetta ‘si’, dunque si è creato un dataset che contiene il 35% di record con etichetta ‘si’, questo ha fatto si che le prestazioni aumentassero in modo significativo.
Gli algoritmi di classificazione utilizzati sono stati: Decision Tree, Naive Bayes, Support Vector Machine e Support Vector Machine con discesa del gradiente, KnearestNeighbor e, per finire, alcuni ensamble learning, quali: RandomForest, Voting e Bagging. I valori di picco per le misure calcolati nella cross validation sono:
• Accuracy: 78% (tutti i classificatori a meno del Decision Tree (76%))
• Precision: 76% (Voting e Knn)
• Recall: 60% (RandomForest e SVC)
• F-Measure: 65% (tutti i classificatori meno DecisionTree e NaiveBabyes)
Per quando riguarda le prestazioni medie, il bagging sembra essere il più equilibrato.
Si è provato ad applicare anche le tecniche di clustering sui dati, dapprima si è provato ad identificare le aree geografiche, ma la mancanza di informazioni non ha permesso una buona clusterizzazione. Si è allora provveduto a provare di modellare i punti del dataset preprocessato in cluster, così da poter applicare al risultato alcuni algoritmi di classificazione nella speranza di poter migliorare gli scores. Putroppo la clusterizzazione non ha dato i risultati sperati. Infine la rete neurale, considerando il dataset bilanciato, presenta più o meno gli stessi valori dei classificatori migliori, però come output fornisce un'informazione aggiuntiva che è la probabilità. Proprio grazie alla probabilità fornita dalla rete che si può modellare ulteriormente il modello, infatti cambiando il risultato della classificazione basandosi sui valori basis che presenta la probabilità, è possibile ottenere dei valori di recall pari al 68%, anche se la presione diminuisce (ma di poco). Nessun classificatore oltre alla rete permette di fare questo. Concludo dicendo che la rete neurale è più versatile e permette all'utente di avere una misura di quanto sia affidabile quella determinata calssificazione.